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INTRODUCTION 


Bak,  Tang,  and  Wiesenfeld  (ref  1)  introduced  the  notion  of  self-organized  critical  (SOC)  systems  to 
provide  a  consistent  explanation  for  the  fractal  spatial  structures,  power  law  distributions,  and  flicker 
noise  commonly  observed  in  spatially-extended  driven  dissipative  dynamical  systems  and  suggested  that 
sand  piles  were  examples  of  such  systems.  They  proposed  that  the  process  of  adding  grains  of  sand  at 
random  locations  within  a  grid  produces  a  sand  pile  whose  geometric  configuration  evolves  toward  a 
critical  form.  The  critical  state  is  an  attractor  for  the  dynamics  of  the  system.  The  implications  of  the 
SOC  concept  for  the  distributions  of  avalanches  in  "sand  pile  numerical  models"  have  been  extensively 
explored  (refs  1-4).  It  was  demonstrated  that  avalanche  sizes  in  finite  sand  pile  models  had  the  form  of 
"cutoff  hyperbolic  distributions"  consistent  with  the  predictions  of  Bak,  Tang,  and  Wiesenfeld  (ref  1). 
These  models  demonstrate  how  complex  fractal  surfaces  evolve  from  simple  system  dynamics.  The  state 
of  the  sand  pile  at  any  time  is  given  by  the  height  of  the  pile  h(x,y)  at  points  {x,y }  on  a  two-dimensional 
lattice.  The  automata  are  governed  by  rules  that  define  the  conditions  for  site  stability  and  the  nature  of 
the  local  rearrangements  that  occur  when  a  site  becomes  unstable. 

Kadanoff  et  al.  (ref  2)  showed  that  the  distribution  of  avalanche  sizes  that  occurs  when  an  SOC  state 
is  perturbed  can  be  fit  over  extensive  ranges  by  power  laws.  Jensen  et  al.  (ref  3)  and  Meisel  and  Cote 
(ref  5)  demonstrate  that  cutoff  power  law  distributions  of  avalanche  sizes  give  rise  to  flicker  noise  in  the 
power  spectral  density.  Thus,  time-scale  invariance  (flicker  noise)  over  increasing  time  scales  and  length- 
scale  invariance  (fractal  scaling)  over  increasing  spatial  scales  develop  as  the  model  is  "driven"  through  a 
sequence  of  subcritical  self-organized  states.  The  range  of  fractal  scaling  correlates  with  the  size-effect 
parameters  in  the  cutoff  power  law  avalanche  size  distributions  and  in  the  critical  state,  fractal  scaling 
may  be  observed  over  essentially  all  scales. 

However,  a  quantitative  analysis  of  the  geometric  scaling  properties  of  the  evolving  automata 
structures  has  not  previously  been  reported.  The  limited  point  sets  generated  in  typical  sand  pile 
simulations  are  insufficient  for  the  application  of  standard  fractal  analysis  techniques.  Therefore,  a  new 
procedure  for  extracting  fractal  parameters  has  been  developed  (ref  6).  This  new  technique  is  inherently 
parallel  and  its  fast  convergence  makes  it  well  suited  to  the  study  of  limited  fractal  subsets.  The 
algorithm  has  been  implemented  on  a  multicomputer  environment  using  a  subset  of  TROLLIUS  (Ohio 
State  University  and  Cornell  Research  Foundation)  called  LAM  (ref  7).  LAM  is  a  programming 
environment  and  development  system  for  a  message  passing  multicomputer  comprised  entirely  of  UNIX 
workstations. 

THE  SAND  PILE  DYNAMICAL  MODELS 

The  cellular  automata  employed  in  this  study  are  based  on  a  two-dimensional  regular  lattice  of  cells 
using  the  rules  similar  to  those  given  for  model  2  of  Reference  2.  The  basic  variable  is  z(ij )  where  (ij) 
are  the  spatial  indices  and  z  represents  the  height  of  the  lattice.  The  system  is  initialized  to  a  planar 
surface,  z(ij)  =  0.  Particles  are  added  at  random  locations  z(ij)  such  that  z(ij)  — >  z(ij)  +  1  as  long  as 
all  the  sites  are  stable.  After  a  particle  is  added  and  after  particle  rearrangements,  the  stability  of  all 
lattice  sites  is  determined. 
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The  stability  of  the  site  [ij)  is  determined  by  z(ij )  and  the  nearest  neighbor  z  values.  In  the  model 
treated  here,  a  site  becomes  unstable  when 


*=i 

£{(z(t,y) - z(i,j-k))H(z(i,j) - z(i,j-k ))  +  (z(i,j) - z(i - k,j))H(z(i,j ) - z(i-k,j))}  > ac 
*=-i 

where  H  is  the  Heaviside  function  and  oc  is  the  stability  parameter. 

If  site  {ij}  becomes  unstable,  z(i,f)  is  reduced  by 


*=i 

^  [H(z(i,  j)  -  z(i,  j-k))  +  H(z(i,j)  -  z(i  - kj ))} 

*=- 1 

and  the  height  of  the  corresponding  nearest  neighbor  cells  is  incremented  by  1  if  and  only  if  their  z-value 
is  less  than  z{i,f).  Note  that  in  Reference  3,  the  stability  of  a  site  is  determined  by  all  nearest  neighbors 
(so  that  a  higher  z  neighbor  site  tends  to  stabilize  a  given  site),  and  that  particles  are  allowed  to  slide 
"uphill."  If  all  sites  are  stable,  a  new  particle  is  added  to  the  system. 

The  boundaries  are  treated  as  cells  with  z(x,y)  =  0.  The  cellular  automata  rule  is  applied  recursively 
to  the  cells  whose  state  is  affected  by  the  unstable  site,  and  the  diffusion  process  continues  until  there  are 
no  more  unstable  sites.  No  new  particles  are  added  to  the  lattice  until  the  lattice  stabilizes.  The  number 
of  sites  that  is  changed  as  the  lattice  reorganizes  after  a  particle  is  added  is  the  size  of  the  resulting 
avalanche,  and  the  total  number  of  particles  that  has  been  added  to  the  system  is  denoted  N. 

D  BY  TRIANGULATION 

The  dynamical  model  under  investigation  yields  a  single- valued  sand  pile  surface  z(i,  j).  It  is, 
therefore,  possible  to  define  a  hierarchy  of  approximations  of  the  surface  in  terms  of  triangles  of  different 
sizes  in  a  simple  unambiguous  way.  The  extent  of  the  lattice  models  studied  has  been  adjusted  to 
facilitate  the  triangular  approximations.  Taking 

L  =  N  -2n  +  1  for  integer  n 

enables  one  to  define  sets  of  triangles  whose  vertices  coincide  with  subsets  of  points  on  the  lattice  and 
whose  projections  in  the  xy-plane  are  right  triangles  of  side  lengths  T,  which  we  refer  to  as  "yardsticks," 
given  by 


Y(m)  =  2m  form  =  {0,1,.. .n] 
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The  small  Y(m)  approximations  to  fractal  surfaces  scale  approximately  according  to 

dln(A(”0)  _ .  2-D 

dln(Y(m)) 

where  the  area  A(m)  is  the  surface  area  for  the  triangulation  based  on  Y(m)  by  Y(m)  cells,  and  D  is  a 
fractal  dimension.  A(m)  can  be  expressed  in  terms  of  individual  triangle  areas 

A(m)  =  XA.  (™) 

Note  that  the  smallest  yardstick  in  the  set  is  Y(0)  =  1,  which  is  as  far  as  one  can  go  in  the  analysis  of  a 
surface  defined  on  a  discrete  lattice  of  points. 

A  detailed  description  of  the  algorithm  has  been  published  (ref  6).  The  reliability  of  the  algorithm  has 
been  successfully  tested  with  Euclidean  and  "Brown"  surfaces  (ref  8).  Brown  constructions  were 
selected  for  the  tests  because  elevations  are  single-valued,  as  required  by  the  algorithm,  and  because  D 
can  be  adjusted  to  have  any  value  between  2  and  3. 

LAM  IMPLEMENTATION 

The  triangulation  algorithm  is  capable  of  reliably  determining  fractal  scaling  parameters  of  relatively 
sparse,  nonhomogeneous  data  sets.  However,  the  generation  of  local  D- maps  is  numerically  expensive. 
Therefore,  a  suite  of  homogeneous  parallel  procedures,  which  exploits  the  inherent  parallelism  of  the 
computational  problem,  has  been  developed  to  map  the  triangulation  algorithm  to  a  multicomputer 
platform  using  a  subset  of  TROLLIUS  (product  of  Ohio  State  University  and  the  Cornell  Research 
Foundation)  called  LAM  (ref  7)  for  interprocess  communication.  The  multicomputer  is  restricted  to  a 
network  of  machines  that  run  LAM  as  a  native  application  under  UNIX.  LAM  generally  supports  only 
general  purpose  UNIX  machines  connected  via  a  local  area  network  or  the  internet.  LAM  is  a  node- 
oriented  computing  environment  that  uses  a  unique  identifier  assigned  to  each  node  (nodeid)  as  the 
primary  synchronization  for  communication.  The  nodes  are  usually  fully  connected  (maximum  1-hop 
distance)  since  the  network  is  generally  a  shared  resource.  In  our  parallel  implementation  of  the 
triangulation  algorithm,  a  master  process  assigns  equal-sized  regions  of  the  structure  to  slave  processing 
nodes  for  computation  of  local  D-map  subsets.  The  processing  nodes  return  the  results  to  the  master 
process,  where  the  D-map  is  assembled  and  the  slave  processes  are  assigned  new  regions  to  analyze. 
Although  this  approach  does  not  exploit  the  communication  advantages  of  a  fully  connected  network,  the 
load  balancing  resulted  in  a  performance  improvement  directly  proportional  to  the  number  of 
workstations.  This  is  likely  due  to  the  high  computation/communication  ratio  and  the  relatively  small 
(15)  number  of  workstations  making  up  the  multicomputer. 
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RESULTS 


Results  were  obtained  for  structures  with  oc  =  20  through  90  in  increments  of  10  on  513  x  513 
lattices  for  extensive  ranges  of  N.  We  refer  to  the  configuration  obtained  for  oc  =  20  and  N  =  7x10** 
particles  as  the  oc  =  20  structure,  and  that  for  oc  =  40  and  N  =  7xl08  as  the  oc  =  40  structure. 

Figure  1  shows  interpolated  shading  representations  of  central  40  x  40  regions  of  the  cc  =  20  and  the 
Gc  =  40  structures. 

Figure  2  plots  D  versus  cc  for  central  regions  of  "close-to-critical"  sand  pile  structures.  As  Gc 
increases,  larger  N  values  are  required  to  reach  critical  sand  pile  configurations.  However,  if  N  is  chosen 
too  large,  boundary  effects  "diffuse"  into  the  region  being  analyzed.  Thus,  the  irregularities  in  the  D 
versus  oc  results  are  consequences  of  varying  degrees  of  convergence  toward  critical  configurations 
rather  than  characteristics  of  the  model. 

Figure  3  is  a  typical  central  cross  section  of  the  oc  =  20  structure  showing  the  boundary-effected 
regions  and  the  relative  variations  in  z.  The  boundary-effected  region  has  diffused  into  the  central  257  x 
257  square;  however,  the  central  129  x  129  square  has  not  yet  been  impinged  upon  by  boundary  effects. 

Figure  4  shows  distributions  of  local  D  values  (obtained  by  analysis  of  the  corresponding  D-maps)  for 
the  oc  =  20  and  40  structures.  Boundary  effects  are  represented  by  the  local  D  values  near  2.0  in  the 
central  257  x  257  region  for  the  oc  =  20  structure.  These  are  consequences  of  the  fact  that  the  steep 
"walls"  seen  in  Figure  3  are  nearly  flat.  The  effect  of  varying  oc  for  regions  free  of  boundary  influences  is 
apparent  in  the  results  presented  for  the  central  129  x  129  regions:  Generally,  as  oc  increases,  the  mean 
value  of  D  and  the  breadth  of  the  distribution  of  local  D  values  increase. 

Figure  5  shows  how  the  structure  with  cc  =  90  passes  through  subcritical  regions  as  more  particles 
are  added.  The  useful  scaling  range  in  the  least  squares  fit  (the  range  over  which  fractal  scaling  is 
obtained)  increases  from  3  doublings  at  N  =  107  to  6  doublings  at  N=  109.  The  range  of  fractal  scaling 
in  the  subcritical  regions  correlates  with  the  magnitude  of  the  cutoff  in  the  size-effect  limited  hyperbolic 
distribution  of  avalanche  sizes. 

SUMMARY 

The  simulations  demonstrate  that  complex  SOC  systems  can  arise  from  the  simple  local  dynamics  of 
cellular  automata.  In  the  present  case,  a  scale-invariant  structure  with  correlating  size-effect  limited 
hyperbolic  distributions  of  avalanches  evolves  from  a  planar  surface  in  a  system  governed  by  simple 
uniform  local  dynamical  rules.  The  diffusion  process  leads  to  uniform  scaling  well  beyond  the  nearest 
neighbor  distances  that  are  directly  involved  when  the  system  is  perturbed. 
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A  new  analytic  technique,  based  on  triangular  tessellations  of  single-valued  surfaces  in  3-space,  has 
been  applied  to  extract  the  fractal  scaling  parameters  from  the  relatively  sparse  data  sets  generated  by  the 
present  simulations.  Since  the  algorithm  is  numerically  expensive  to  apply  to  the  determination  of  local 
D-maps  of  heterogeneous  surfaces,  the  computational  problem  has  been  mapped  to  a  multicomputer 
platform  comprised  of  UNIX  workstations  using  a  subset  of  TROLLRJS  called  LAM  for  interprocess 
communication.  Computation  of  large  sets  of  local  D-maps  is  presently  impractical  without  the  use  of  a 
highly  dedicated  computer  ensemble. 
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Figure  la.  Interpolated  Shading  Representation  of  Central 
Region  of  Sand  Pile  Structure  for 

Gc  =  20  and  N  =  7x10  8. 


Figure  lb.  Interpolated  Shading  Representation  of 
Central  Region  of  Sand  Pile  Structure  for 
Oc  =  40  and  N  =  7x10  8. 


oc  =  20  and  N  =  7x10 
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Figure  4.  Distributions  of  Local  D  Values. 


Figure  5.  Area  vs.  Yardstick  plots  for  oc  =  90.  The  lines  are 
least  squares  fit  to  the  first  4  points  for  the  N=  10^  structure 
and  to  the  first  7  points  for  the  N=  10^  structure. 
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AMSTA-AR-AES,  BLDG.  321  1 

AMSTA-AR-AET-O,  BLDG.  183  1 

AMSTA-AR-FSA,  BLDG.  354  1 

AMSTA-AR-FSM-E  1 

AMSTA-AR-FSS-D,  BLDG.  94  1 

AMSTA-AR-IMC,  BLDG.  59  2 

PICATINNY  ARSENAL,  NJ  07806-5000 

DIRECTOR 

U.S.  ARMY  RESEARCH  LABORATORY 
ATTN:  AMSRL-DD-T,  BLDG.  305  1 

ABERDEEN  PROVING  GROUND,  MD 
21005-5066 


DIRECTOR 

U.S.  ARMY  RESEARCH  LABORATORY 
ATTN:  AMSRL-WT-PD  (DR.  B.  BURNS)  1 

ABERDEEN  PROVING  GROUND,  MD 
21005-5066 

DIRECTOR 

U.S.  MATERIEL  SYSTEMS  ANALYSIS  ACTV 
ATTN:  AMXSY-MP  1 

ABERDEEN  PROVING  GROUND,  MD 
21005-5071 


COMMANDER 

ROCK  ISLAND  ARSENAL 

ATTN:  SMCRI-SEM  1 

ROCK  ISLAND,  IL  61299-5001 

MIAC/CINDAS 
PURDUE  UNIVERSITY 

2595  YEAGER  ROAD  1 

WEST  LAFAYETTE,  IN  47906-1398 

COMMANDER 

U.S.  ARMY  TANK-AUTMV  R&D  COMMAND 
ATTN:  AMSTA-DDL  (TECH  LIBRARY)  1 

WARREN,  MI  48397-5000 

COMMANDER 

U.S.  MILITARY  ACADEMY 

ATTN:  DEPARTMENT  OF  MECHANICS  1 

WEST  POINT,  NY  10966-1792 

U.S.  ARMY  MISSILE  COMMAND 
REDSTONE  SCIENTIFIC  INFO  CENTER  2 

ATTN:  AMSMI-RD-CS-R/DOCUMENTS 
BLDG.  4484 

REDSTONE  ARSENAL,  AL  35898-5241 
COMMANDER 

U.S.  ARMY  FOREIGN  SCI  &  TECH  CENTER 
ATTN:  DRXST-SD  1 

220  7TH  STREET,  N.E. 

CHARLOTTESVILLE,  VA  22901 

COMMANDER 

U.S.  ARMY  LABCOM,  ISA 

ATTN:  SLCIS-IM-TL  1 

2800  POWER  MILL  ROAD 

ADELPHI,  MD  20783-1145 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH,  DEVELOPMENT,  AND  ENGINEERING  CENTER, 
BENET  LABORATORIES,  CCAC,  U.S.  ARMY  TANK-AUTOMOTIVE  AND  ARMAMENTS  COMMAND, 

AMSTA-AR-CCB-O,  WATERVUET,  NY  12189-4050  OF  ADDRESS  CHANGES.  _ 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST  (CONT’D) 


NO.  OF 
COPIES 


COMMANDER  WRIGHT  LABORATORY 

U.S.  ARMY  RESEARCH  OFFICE  ARMAMENT  DIRECTORATE 

ATTN:  CHIEF,  IPO  1  ATTN:  WL/MNM 

P.O.  BOX  1221 1  EGLIN  AFB,  FL  32542-6810 

RESEARCH  TRIANGLE  PARK,  NC  27709-2211 

WRIGHT  LABORATORY 

DIRECTOR  ARMAMENT  DIRECTORATE 

U.S.  NAVAL  RESEARCH  LABORATORY  ATTN:  WL/MNMF 

ATTN:  MATERIALS  SCI  &  TECH  DIV  1  EGLIN  AFB,  FL  32542-6810 

WASHINGTON,  D.C.  20375 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH,  DEVELOPMENT,  AND  ENGINEERING  CENTER, 
BENfiT  LABORATORIES,  CCAC,  U.S.  ARMY  TANK -AUTOMOTIVE  AND  ARMAMENTS  COMMAND, 
AMSTA-AR-CCB-O,  WATERVLIET,  NY  12189-4050  OF  ADDRESS  CHANGES. 


NO.  OF 
COPIES 


1 


1 


